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Abstract We hereby study the stability of a massless probe orbiting around an oblate 
central body (planet or planetary satellite) perturbed by a third body, assumed to lie 
in the equatorial plane (Sun or Jupiter for example) using an Hamiltonian formalism. 

We are able to determine, in the parameters space, the location of the frozen orbits, 
namely orbits whose orbital elements remain constant on average, to characterize their 
stability/unstability and to compute the periods of the equilibria. 

The proposed theory is general enough, to be applied to a wide range of probes 
around planet or natural planetary satellites. 

The BepiColombo mission is used to motivate our analysis and to provide specific 
numerical data to check our analytical results. 

Finally, we also bring to the light that the coefficient J2 is able to protect against 
the increasing of the eccentricity due to the Kozai-Lidov effect. 

Keywords Methods; analytical study • Stability ■ Long-term evolution ■ Kozai 
resonances ■ Frozen Orbit equilibria 



1 Introduction 

BepiColombo (MPO and MMO orbiters) is a joint European and Japanese space agen- 
cies space mission aimed at studying the planet Mercury. The MPO (Mercury Planetary 
Orbiter) will be brought into a polar elliptical orbit around Mercury with an inclina- 
tion of 88-90°, an eccentricity of 0.1632 and a semi-major axis of 3 394 km. The MMO 
(Mercury Magnetospheric Orbiter) will also be brought into a polar elliptical orbit with 
an eccentricity of 0.6679 and a semi-major axis of 8 552 km. 

Actually polar orbits are very interesting for scientific missions to planetary satel- 
lites (with near polar low-altitude) or to planet (with high-eccentric high-altitude). The 
orbital dynamics of such space probes is governed by the oblateness ( J2 effect) of the 
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central body around which the space probe is orbiting and the gravity field from the 
third body. A well-known effect of the third-body perturbation is the change in the 
stability of circular orbits related to orbit inclination. This effect is a natural conse- 
quence of the Kozai-Lidov resonance (Kozai 1962; Lidov 1963). The final fate of such 
a satellite is the collision with the central body. Therefore the control of the orbital 
eccentricity leads to the control of the satellite lifetime. 

Scheeres ct al (2001) studied near-circular orbits in a model that included both the 
third body's gravity and J2. In addition San- Juan et al (2006) studied orbit dynamics 
about oblate planetary satellites using a rigorous averaging method. Paskowitz and 
Scheeres (2006) added the effect of the coefficient J3. These authors mainly focused 
their attention to an orbiter around planetary satellites especially for Europa orbiter. 
So they did not take into account the eccentricity of the third body and they detailed 
the near-circular orbits. 

Our purpose is to build a simplified Hamiltonian model, as simple as possible, which 
will reproduce the motion of probes orbiting an oblate central body also taking into 
account the third body effect. Especially we are looking for the conditions that give 
rise to frozen orbits. Frozen orbits are orbits that have orbital elements constant on 
average. These particular orbits are able to keep constant the eccentricity. Therefore in 
a neighbourhood of these orbits there is a stability area where even a limited control 
couki be used to avoid the crash onto the central body. 

Beside the oblateness of the central body and the gravity effect of the third body, 
our averaged model takes into account also the eccentricity of the orbit of the third 
body. Moreover let us observe that our results arc given in closed form with respect to 
eccentricity and inclination of the probe, namely we do not perform any power series 
expansion; therefore, our theory applies for arbitrary eccentricities and inclinations of 
the space probe, and is not limited to almost-circular orbits. We can thus conclude 
that the theory is general enough to be applied to a wide range of probes around a 
planet or around a natural planetary satellite and, can be formulated and presented in 
a general way that allows extension of the results to other cases. 

The Mercury orbiter mission (BepiColombo) is used to motivate our analysis and 
to provide specific numerical data to check our analytical results. 

We are able to provide the location of frozen orbits and study their stability as a 
function of the involved parameters, using implicit equations and graphics. Finally we 
give the analytical expressions of the periods at the stable equilibria. 

The analytical results are verified and confirmed using dedicated numerical simu- 
lations of the whole model. 

To conclude, we discuss the effect of the protection of J2 on the increase of the 
eccentricity due to Kozai-Lidov effect and the apparition of an asymmetry caused by 
the addition of the coefficient J3. 

2 Motivation: numerical exploration 

For the purpose of our study, we consider the modeling of a space probe subjected 
to the influence of Mercury's gravity field (in the following sections Mercury will be 
denoted by "central body") and the gravitational perturbations of the Sun (noted 
"third body" ) as well as to the direct solar radiation pressure without shadowing effect. 
As a consequence the differential system of equations is given by 



r = Vpot + rQ+ rrp , 



(1) 
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where i-pot is the acceleration induced by Mercury's gravity field, vq is the acceleration 
resulting from the gravity interaction with the Sun and rrp is the acceleration due to 
the direct solar radiation pressure. 

It is worth noting that we modelise the gravity potential of central body only using 
the J2, C22 and J3 coefficients. In our implementation, we choose the high accurate 
Solar System ephemeris given by the Jet Propulsion Laboratory (JPL) to provide the 
positions of the Sun (Standish 1998). We adopt the variable step size Bulirsch-Stoer 
algorithm (sec e.g. Stoer and Bulirsch 1980) to numerically integrate the differential 
equation (1). L(;t us note that, for the purpose of validation, we also use a second 
numerical integrator D0P853 (an explicit Runge-Kutta method of order 8(5,3) with 
stepsize control due to Dormand & Prince (Hairer et al 1993)). 

In Figure 1 we report the results of a numerical integration of the system of equa- 
tions (1) for a set of 19 600 orbits, propagated over a 200 years time span with a 
entry-level step size of 300 seconds. We consider a set of initial conditions defined 
by an eccentricity grid of 0.005 and a semi-major axis grid of 35 km, spanning the 
[2600, 7600] km range. The other fixed initial conditions are iq = 90° for the inclina- 
tion, Qq = 67.7° ujQ = —2° for the longitude of the ascending node and the argument of 
periherm, respectively; Mq = 36.4° the mean anomaly at epoch fixed at 14 September 
2019. The area-to-mass ratio A/m = O.Olm^/kg. These values have been fixed by the 
initial conditions of BepiColombo mission found in Garcia et al (2007). 

We show the amplitude of the eccentricity (that is the difference between the max- 
imum and minimum eccentricity reached during the integration) of each orbit in the 
left panel of Figure 1. For each orbit, we also calculate using the Numerical Analysis 
of Fundamental Frequencies, for short NAFF (Laskar 1988, 2005), the fundamental fre- 
quency (noted v) of the evolution of the eccentricity vector (ecos w, esinw). We plot 
the logarithm of the second derivative (noted \og{55v)) of this frequency in the right 
panel of Figure 1, namely an indicator of the diffusion in the frequency space, hence the 
regularity of the orbit. For more details concerning this use of frequency analysis, see 
Lemaitre et al (2009) where the frequency analysis has been used to study resonances 
in Geostationary Earth Orbits. 

First, let us observe that the white zone in Figure 1 corresponds to orbits that 
crash onto central body's surface. Second we distinguish a curve where the variation 
of the eccentricity amplitude is null (dashed black line) . On the second derivative plot 
(right panel) we also distinguish on the left of the dashed black line (null-variation 
of eccentricity) a larger value of the log of the derivative that could correspond to a 
separatrix. 

These structures will be analyzed and explained using a simplified model, that 
takes into account the central body attraction with the J2 harmonic coefficient and 
the third body gravitational effect. We observed that the solar radiation pressure does 
not play any role in these structures, hence this effect will be absent in the simplified 
model. 



3 The Hamiltonian Formalism 

The aim of this section is to introduce the Hamiltonian (2) already found in Tremaine 
et al (2009). Kepler's Hamiltonian describing the motion of a test particule orbiting an 
isolated point mass M is 

1 2 GM GM 
^K = i:v = 

2 r a 
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Fig. 1 The eccentricity computed as a function of the initial eccentricity eo and the initial 
semi-major axis ao- The equations of motion include the central body attraction, the harmonics 
J2, C22, Js, the solar interaction as well as the perturbing effects of the solar radiation pressure 
{A/m = O.Olm^/kg). The eccentricity step is 0.005 and the semi-major axis step is 35 km. 
The initial conditions are io = 90°,f2o = 67.7°, a;o = -2° and Mq = 36.4°. The integration 
time is 200 years from epoch fixed at f4 September 20f9. The patterns have been obtained 
by plotting the amplitude of variation of the eccentricity and log(5(5i/) respectively in left and 
right panel. 



where G is the gravitational constant, r is the planetocentric position of the particule, 
V — r, r = \r\ and a is the semi-major axis of the particule. 

One can introduce the quadrupole potential arising from an oblate planet ( "central 
body") that is 



GMJ2Rp 
2r^ 



|^3(r ■ rip 



where rip is the unit vector oriented to central body's spin axis (see Figure 2). M, 
Rp and J2 are respectively the mass, the radius and the oblateness coefficient of the 
central body (planet or natural satellite). 

We assume that r <^ a^^ (where the subscript 3^ is related to the third body) and 
we average over the third body orbital period. So, we obtain the quadrupole in terms 
of the third body gravitational effect 



4a3,(l-e2j3/2 



where n^^ is the normal to the central body orbit. M^^ 



and are respectively the 



mass, the semi-major axis and the eccentricity of the third body. This quadrupole term 
takes into account the eccentricity (63^) of the third body (e.g. Sun or Jupiter) around 
the central body (planet or natural satellite). Let us stress the fact that Paskowitz 
and Scheeres (2004); San- Juan et al (2006); Scheeres et al (2001) do not include this 
eccentricity factor in their formulation, while for a Sun-Mercury-orbiter application, 
this will be an important contribution. 

We then average over the Keplerian orbit of the test particule described by the 
following elements: a semi-major axis a, an eccentricity e, and an orientation specified 
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by the unit vectors n along the angular momentum vector, u toward the pericenter 
and V = n X u. We have (Brouwer and Clemence 1961) 



<{vuf> = am+2e'), < {r ■ vf > = a^'^ - '^e^) , 



(r ■ uf \ _ / (r ■ vf 



2a3(l-e2)3/2- 

where < > denotes the average over M, the mean anomaly of the orbit. 

r ■ A ^ [gM J2RI , Mg.a^ 

Let J = V 1 - n , e = e-u, r = 1 / ^^f, e,^ = — — and e.,, = 5— ^ryr 

where e is the eccentricity vector and Ej^ >Q, £31, > 0. We finally define a dimensionless 
(divided by GM/a) Hamiltonian 

^' = -\ + _'e2)5/2 [1 - - 3(J • np)'] + [5(e -n^.f- (j ■ n,,f - , (2) 

That describes the secular equations of motion of a test particule around an oblate 
central body perturbed by the third body gravitational effect. Let us summarize the 

used assumptions: 

1. the precession rate of the central body spin due to third body gravity is negligible; 

2. the satellite is a massless test particule; 

3. the third body is far enough from the central body such that the third body gravity 
can be approximated by a quadrupole; 

4. the satellite is far enough from the central body so that the potential from the 
central body can be approximated as a monopole plus a quadrupole; 

5. the perturbing forces {$,2 + ) are weak enough so that the secular equations of 
motion can be used to describe the orbital motion; 

6. there are not resonant relations in mean motions between the frequencies of the 
satellite and the frequencies of the central body. 

Let us remark that San-Juan et al (2006) already studied the orbit dynamics about 
planetary satellites using an extensive averaging method based on the Lie transforms to 
obtain averaged equations involving higher orders whose result is the introduction of an 
asymmetry for direct and retrograde satellite. Our simplified model will not be able to 
capture this asymmetry because the resulting Hamiltonian (3) will be symmetric in the 
satellite inclination; thus direct and retrograde satellites will have the same behavior. 

Let us now make some assumptions suitable in the case of a non-inclined central 
body orbit (e.g. Sun-Mercury-orbiter system or Jupiter-Europa-orbiter system). We 
hereby consider an equatorial third body, thus rip = ng^. We also set G — vl— 
and H = G cosi where j ■ rip — Vl — e2 cosi. To eliminate an extra parameter, we 
divide the Hamiltonian by the coefficient Sj^ and we introduce the coefficient 7 



ej. Ma3,(l-e2j3/2 j^i?: 



2 • 



In Figure 2, we represent the geometry for the general problem (on the left) and for 
our simplified one (on the right). 
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The averaged Hamiltonian is then 



K.' 



1 

4G3 



3^ 



+ 



+ 



. 2 

sm u! - 



i/^ - 2 + 2G^ 



37 
8 



5(1 - 



1 - 



El 

G2 



sin^ — H'^ — 2 



2G^ 



(3) 

This Hamiltonian is independent of the ascending node Q. If we take 7 = (£3^ = 0, 
namely we take into account only the oblateness effect), we have the well-known circular 
dynamics of the eccentricity vector due to the J2 coefficient with an elliptical fixed 
point in the semi-equinoctial elements {k,h) = (Vl — G2 cos a;, Vl — G2 sinctj). If we 
take 7 — >■ 00 {e,2 ~ i.e. only the third body contribution does matter), we find 
the Kozai-Lidov Hamiltonian which we find in a similar formalism in Paskowitz and 
Scheeres (2004) (with e^^ — 0). The Hamiltonian (3) (with e^^ — 0) can also be found 
in Scheeres et al (2001). 
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Fig. 3 Relation between 7 and the semi-major axis of a test particule orbiting the central 
body (terrestrial planets or Europa). The third body are respectively the Sun and Jupiter. 
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Table 1 Connection between semi-major axis (km) of the probe around the central body and 
the coefficient 7. The rows "Min." , "Missions" and "Hill" give the values of 7 with respect to 
the semi-major axis respectively equal to the radius of the central body, to one space mission 
and to the radius of the Hill's sphere. The "--" symbol indicates that the semi-major axis is 
lower than the radius of the central body or greater than the radius of the Hill's sphere. 





Min. 


Missions 


Hill 


Particular values 


Mercury 


a 


(km) 


2 439.990 


Messenger 
10136.2 


175 295 


4 350 


5 577 






7 


0.008 


9.9136 


1.533 X lO'' 


1/7 


0.5 




a 


(km) 


6 051.8 


Venus Express 
39176.8 


1 004 270 


9 350 


12 010 






7 


0.008 


184.485 


2.042 X lO'J 


1/7 


0.5 


Earth 


a 


(km) 


6 378.137 


Meteosat 
42164.14 


1 471 506 


36 350 


46 670 






7 


2.38 X 10-5 


0.30107 


1.559 X 10^ 


1/7 


0.5 


Mars 


a 


(km) 


3 396.190 


Mars Express 
9311.95 


982 748 


26150 


33 580 






7 


5.288 X 10-6 


8.195 X 10-4 


1.073 X lO'^ 


1/7 


0.5 


Europa 


a 


(km) 


1565.0 


EJSM/JEO 
3 222 


13 529 










7 


1.153 


42.646 


5.568 X 10-* 


1/7 


0.5 



For illustration, we respectively show in the Table 1 and draw in Figure 3 the value 
of the coefficient 7 with respect to the semi-major axis for a probe around a terrestrial 
planet and around Europa. This coefficient can be related to other parameters used 
in the literature. For example, it can be linked to the coefficient /3 in San- Juan et al 
(2006) or to the coefficient e used in Scheeres et al (2001). 



4 Secular Equations of Motion 

From the doubly averaged Hamiltonian (3), we obtain the equations of motion: 



H = 



37 



5 I ^ - G ] sin^ LU + 2G 



G = ( 1 — G ) 1 — T7 sm bj cos uj 

4, \ G2 ' 



3 

+ 4G^ Tg^^^ 



(4) 

(5) 
(6) 



Developing these equations in eccentricity up to second order, we can obtain the equa- 
tions of Scheeres et al (2001). In the following, we will adopt a complementary approach, 
keeping functions of eccentricity and inclination, without power series developments, 
in such a way that our results hold for any arbitrary eccentricities and inclinations. 
From the previous set of equations, we observe that 77 is a constant of motion, and 
cos^ I too, as in the Kozai-Lidov effect (Kozai 1962; Lidov 1963). Besides, let 
us remark that < G < 1, thus < 77 < G < 1 and moreover 7 > 0. The first equation 
(4) is equal to zero only for 77 = corresponding to exact polar inclination. Moreover 
the ascending node does not affect any of the other orbital elements. The last equation 
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(6) equals to zero for G = 1, u = kTi/2, fc G N or _ff = G, namely i = 0°, that is the 
planar case. The third equation (5) could equal to zero for oj = 0, tt or a; = ±7r/2. We 
analyze these equations in next section to find the equilibria. 



5 Frozen Orbit Solutions 

A frozen orbit is characterized by no secular change in orbital eccentricity and argument 
of pericenter. It has constant values of e, i and w on average, this results in fixed 
geometrical size and locations, apart from short period oscillations. 

We already observed that equilibria appear when G = lora' = 0,7rora) = ±7r/2. 
We separatly deal with these three different cases. For each of them, we give the number 
of equilibria, the conditions of existence and we calculate their stability. 

Wc do not deal with the singularity G = Q (<4> e = 1) because it corresponds to 
an escape of the orbiter. We will show that the equilibrium G = 1 always exists. So to 
begin, we deal with the non-circular case G =^ 1 (eccentricity =^ 0). 



5.1 Non-circular case G ^ 1 (eccentricity ^ 0) 

5.1.1 Vertical equilibria - Kozai-Lidov equilibria: cosw = w = ±7r/2 
The conditions to simultaneously equal to zero the equations (5) and (6) is: 

f „2_ 1 + 3GS 

^ - 5 1 + G37 (7) 

I cos w = 

Because < G < 1 then this equation implies that 

< . (8) 

Lot us observe that this is also the value for which one real root does exist. If this 

condition is violated then no real root exists. 

Actually we determine a region given by the implicit equation 

864 0001/^^7*^ + (^2 963 520//^^ - lQ2AH^°^'y'^ 

(^1 512 630 il* - 13 9651/^ - 22 235 661 j 7^ + 12 = Q 
and H"^ <^ 

where it is possible to find three real roots. We will show that these three reals roots 
appear for eccentricities larger than 0.996 59. Being a case close to an escape of the 
orbiter, we will leave to section 7.3.3 a discussion of this "local deformation". 

If the oblateness term is neglected (sja ~ <^ 7 — > 00), the existence condi- 
tion becomes independent of the physical parameter and reduces to sin^ z < | or 

arccos ~ 39.23° < i < 144.77° which corresponds to Kozai-Lidov critical inclina- 
tion. 
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Wc also analyze the stability of these equilibria (7). The Jacobian of the Hamilto- 
nian (3) evaluated at the equilibrium (7) (noted by |Eq.(7) or Gi~i being the value of G 
at the Kozai-Lidov equilibrium) is given by: 



Eq.(7) 



Eq.(7) 



dGduj 



Eq.(7) 



3 
2G5 



15 



G2 



4 V^^G^ 



-jG^ 



kl 



Eq.(7) 

2 + 7G|,-2l7G|i- 127^(31;) 



15 



-^7(1 -G^) 
-^7(1 -Gi) 



£1 
" G2 



Eq.(7) 



(9) 



dojdG 



1 + 
0. 



Eq.(7) 



In the equations (9), the term ^ ^ is always strictly negative (if G < 1). Then 

the equilibrium is a stable point if 



1 

G5 



15 



■El 

G2 



37 
2 



1 + 5 



■El 

G4 



< 



■2l7G|;+157Gfei/f2<0. 



Eq.(7) 



(10) 

This equation (10) is always satisfied for all 7 > 0, < and G^i < 1 {ski > 0). 

Therefore, for these conditions, we have two opposite stable points at a) = ±7r/2 and 
G such that H'^ = ^ TTU^- 

For the inclination t = 90° {H^ = 0) the equilibrium exists for the particular value 
of G = (e = 1). This case is only theoretical and should not be considered, because 
it would correspond to an escape of the orbiter. In Figure 4 we give the location of the 
equilibria G (7) in the parameter space (7,i?^). 



5.1.2 Horizontal equilibri,a: sinw = <^ w = 0, tt 

The conditions to simultaneously equal to zero the equations (5) and (6) are: 



= ^ (1 - 2GS) 
sin a; = 0. 



(11) 



Using "Le theoreme d'algebre de Sturm" (Sturm 1835) (for more explanation see the 
Appendix) we calculate the number of roots (G) in the range < G < 1 of the equation 
(11) as a function of the parameters 7 and . For 7 > this equation has 

• one real root, equal to if = and 7 < 1/2. 

• three real roots (one equal to and the other two opposite) if = and 7 > 1 /2. 

• three real roots (one equal to 1 and the other two opposite) if < < ; 

• five real roots (one equal to 1 and the other ones opposite two by two) if 7 > 1/7 



1^2 < 



(77) 



-2/5 



and — g — ^ ^ Y- 

one real root equal to 1 otherwise. 
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Fig. 4 Values of G at Kozai-Lidov stable (Eqs. 7 and 8) equilibria (vertical equilibria: uj = 
it7r/2) computed as a function of and 7. These equilibria are always stable. The color code 
indicates the value of G at the equilibrium. 



In Figure 5, we give the location of the equilibria G (11) in the space {^,H^). The 
particular case G = 1 will be treated in the next section. We can also analyze the 
stability of these equilibria (11). The Jacobian of the Hamiltonian (3) evaluated at 
the equihbrium (11) (noted by |Eq.(ii) or Ghor, being G^or the value of G at the 
equilibrium) is given by: 



( 



Eq.(ll) 



2G5 V G2 j + 2 



Eq.(ll) 



2Gi 



-l + 7jG 



hor 



Eq.(ll) 



15 L H' 

^7(1-G)(1-^ 



(12) 



Eq.(ll) 



dG du 



^7(l-GLr)(2 + GLr7) 



Eq.(ll) 



du dG 



In the equations (12), the term 
the equilibrium is a stable point if 

1 



lEq.(ll) 



7 > 



G5 



= 0. 

Eq.(ll) 

is always strictly positive (if G < 1). Then 



^5-1 

77 



Using this equation at the equilibrium (11), we obtain conditions for stability of the 
stable point (G / 1 <4> e / 0) 

2/5 



1-27 ,,2 1/1 



and - < 7 . 

7 - ' 



7 V77 



(13) 
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So the condition to have an unstable equihbrium is given by 

1 „2 1-27 „2 1/1 

7<- or H'<^ or > ^ [^) ■ (14) 
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Stable equilibrium 



Unslable equilibrium 
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0.05 
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Fig. 5 Value of G at stable (on the left) (Eqs. 12 and condition 13) equilibrium and unstable 
(on the right) (Eqs. 12 and condition 14) equilibrium, computed as a function of and 7. 
The color code indicates the value of G at the equilibrium (horizontal equilibria: ui = 0, tt). 
The inset on the left panel shows the same plot of left panel but using a wider color code. 



In the Figure 5, we notice that when both unstable and stable equilibrium exist, 
the unstable equilibrium always appears for a value of G lower than that of the stable 
point (i.e. for a value of e greater than the one for the stable point). 



5.2 Circular case G = 1 (eccentricity e = 0) 

For the case G = 1, we can use a canonical transformation to cartesian coordinates 

X = a/ 1 — sin bj y = \/ 1 — G^ cos cj (15) 

The new Hamiltonian is therefore 



K. = 



4 V(l-2;2 -y2)3/2 (1 _ 3.2 _ 2^2)5/5 



37 



5a;^ ( 1 



1 ~ x-' ~ y" 



H"^ - 2x^ ~ 2y^ 



(16) 



for which it is obvious that (0, 0) is always an equilibrium point whose stability can be 
studied computing the second derivatives and evaluate them at this equilibrium: 



12 



dx \x=0=y 4 4 



— I - -(1 

Qy^ \x—Q—y 4 



Si/-" 



37 
2 

= 0. 



. dx dy\x=Q=y dy dx\x=o=y 
So, the condition to have a stabihty point at x = = j/ is 



27 



or 



^2^1 + 37 



5 57 + 5 ' 

and thus the condition to have an unstable point at a; = = y is 

1 ^ ^2 ^ 1 + 37 



57 + 5 



(17) 



(18) 



5.3 Summary of the phase space 

In this section we summarize the various possible phase spaces topologies as a function 
of the parameters. We draw (Fig. 6) the bifurcation lines (conditions 8, 13 and 17) in 
the parameter space (7,//'^). This bifurcation diagram is equivalent to the upper part 
of the bifurcation diagram in San- Juan et al (2006) but here we draw the bifurcation 
lines in the general (not linked to a particular central body) space 

(7,^2). The H'^ = 

(77)~^/'^/7 line stops at the limit 7 = 1/7. For this value, this curve coincides with 
the = (1 — 27)75 condition. For the Jupiter-Europa-orbiter system, the minimum 
value of 7 is 1.153 (Tab. 1). Therefore the phase spaces (A) and (E') do not exist. 



0.6 




Fig. 6 Bifurcation lines and regions in the parameter space {■y,H'^). In the regions (B) and 
(C) we have the same number and the same stability of the equilibria but the phase space 
is topologically different. These regions are separated by the dashed line implicitely given by 
equation (19). 



The region (E') and (F) in magenta color correspond to exact polar orbits (i = 90° 
thus H"^ ^0). 
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For each region, we attribute a letter and we draw (Fig. 7) a generic contour plot 
of the Hamiltonian (3) in the (fc, h, t) physical space. We recall that the motion of 
the inclination i is given by the conservation of the first integral H = G cost. We 
also draw the projection of these phase spaces in the semi-equinoctial elements space 
{k, h) = (Vl — cos w, Vl — sin w). In this phase space, it is easier to bring to the 
fore the stable (green point) and unstable (red cross) equilibria. The (E') phase space 
is trivial, containing only concentric circle in the i = 90° plane, so we do not repoduce 
it. 

In the Figure (7), we notice that the maximum inclination is always reached at 
e = 0. This is explained by the relation = Vl — cos i. This last relation also 
gives a maximum bound onto the eccentricity: e < Vl — H'^. Therefore there are some 
values of H for which the phase space is visibly restricted in eccentricity. Beyond this 
eccentricity, the motion is physically impossible. 

Let us observe that the region near the stable equilibria allows to control the 
variation of the eccentricity even for high eccentricity. We also remark that there are 
"dangerous" portions of phase space such as the region around the 7 = = 1/7 or 
near of the (B)-(C) transition. In these regions the dynamics (in a full model) could 
change strongly for a small variation of (//^,7) or (e,oj). 

The transition between (B) and (C) phase spaces arises when the cncrgj' of the 
separatrix at the (0, 0) equilibrium is equal to the energy of the unstable exterior 
horizontal equilibrium. This condition gives a new "fictitious" bifurcation line (dashed 
line in Figure 6) in the parameter space (7,H^). To find this line, we evaluate the 
Hamiltonian (3) at the unstable equilibrium = ^|-(1 — 2G^7) (Eq.ll and condition 
14) and we denote this value by ICi. Afterward, we evaluate the Hamiltonian (16) at 
the unstable equilibrium (0, 0) (condition 18) and denote the result by JC2 ■ We now 
assume these two equilibria have the same value of Hamiltonian K, and of H^. Then 
we can replace H"^ by ^{1 - 2G^7) in K.2 and we impose the equality between ICi 
and /C2. Therefore we obtain the condition 

nor ' nor nor 

q2 

where G^g^ is the unstable horizontal equilibrium i.e. H = — !i2L(^i _ 2G/jo^7). We 

5 

plot this implicit condition (19) in Figure 6 with a dashed black line. This line joins 
the "(1 - 27)/5" and "(77)-2/^/7" lines at the (7 = 1/7, H"^ = 1/7) point. 

For the particular case 7 — >^ (J2 effect only), we obtain, for all H^, a phase 
space with circular motion of the eccentricity. We see that near to the value = 1/2 
(corresponding to the Molniya^ critical inclination equal to 63.43° with G = 1), the 
phase spaces (A), (D) and (E) always exist until 7 becomes exactly equal to 0. In the 
opposite case, 7 — ^ 00 (third body effect only), the curve = (77)~^/^/7 converges 
to and the curve = {1 + 37)/(57 + 5) converges to 3/5 (corresponding to the 
Kozai-Lidov critical inclination equal to 39.23° with G = 1). Then only the following 
phase spaces are reaUzable: (E) (for H'^ > 3/5), (D) (for < ^ < 3/5) and (F) (for 
= 0) with (F) that degenerates to an unstable point at the center. These three 
phase spaces will be shown in Figure 12. 



^ At this inclination, due to J2 effect, the argument of perigee remains nearly constant for 
a long period of time. Moliiiya orbits are named after a series of Soviet/Russian Molniya 
communications satellites which have been using this type of orbit since the mid 1960s. 
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Fig. 7 Examples of some generic contour plots of the Hamiltonian (3) in (k, h, i) space for 
each region of Fig. 6. The inclination i is given in degrees and the semi-equinoctial elements 
(k, h) are given by (k, h) = — G'' cos uj, Vl — G"^ sintj). The green point and red cross are 
respectively the stable and unstable points. In the polar projection, the radius is the eccentricity 
e and the angle is the pericenter lo in degrees. 
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In Figure 8, we show how the stable and unstable equilibria evolve, appear and 
disappear in each region and during the transition between the regions. We take a 
vertical section in the Figure 6 at 7 = 0.4. This section crosses the regions (A), (C), 
(B), (D) and (E). We draw the value of the eccentricity at the stable (solid color 
lines) and unstable (dashed color lines) equilibria with respect to . The vertical 
dashed black lines mark the boundary of the regions. The numbers give the number of 
equilibria with this value of e. For example, 2 in magenta dashed line means that there 
are two unstable equilibria with the same value of e, respectively for cj = and u = n. 

At the transition between (A) and (C), the central (e = 0) stable point bifurcates 
in two horizontal stable points (e 7^ and cj = 0, tt) and one unstable point (e = 0). 
At the transition between (B) and (D), the two unstable and the two stable horizontal 
(oj = 0, tt) equilibria converge to the same value of e and cancel out. At the transition 
between (D) and (E), the two stable Kozai-Lidov equilibria (lj — ±vr/2) come close 
to and cancel out with the central unstable equilibrium to give one central stable 
equilibrium. We remark that the transition between (C) and (B) is not characterized 
by a change of the equilibria. 



5.4 Period at the equilibrium 



We are now interested in the period of the eccentricity vector at the equilibrium. This 
will be done by linearizing in a neighborhood of the equilibrium. Then the Hamilto- 
nian close to the equilibrium is given by (the subscript eg. means "evaluated at the 
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equilibrium" 

/C = K.eq. + ^1 {G-Geq.) + TT" I (w-Weg.) 

' • leg. 





Geg.)(w-Weg.) + ^ 

This is an harmonic oscillator that can be expressed in action-angle variables {ij), J) 
defined as (at a stable equilibrium, we have ab> 0): 

X=t-V2JcostP and Y=(^V2Jsmtp. 
\ a \ 

Then the frequency at the equilibrium is given by 



Using the equation (20), the periods (t) at the stable equilibria are given by: 
— for horizontal equilibria: G such as Equation (11) and condition of stability (13) 

47r 



V " ^^^) + (1 - G2) (l - ^) 

for vertical (Kozai-Lidov) equilibria: G such as Equation (7) and condition of sta- 
bility (8) 

5 r£zi^_i5^) _!^fi + 5Hj)l,^^(l_G2)fi_ii^ 



— for central equilibrium (e = 0): G = 1 with condition of stability (17) 

^ Stt 

We remind that 7 = e.-ib/e./g a-nd that the equations are dimensionless. Then the periods 

at the equilibria are given by Teq. — \J req. . 

For example, we apply these formula to a Mercury orbiter. The values for Mercury 
are a^^ = 57909176.0 km, e^^ = 0.205 630 69, J2 = 6.0 x 10"^ (Anderson et al 1987) 
and Rp = 2 439.99 km. In the Figure 9 we plot the periods at the equilibria respectively 
for the three cases: 

— on the left panel, the periods at the stable equilibrium with respect to the value 
of 7 and H^. The color code indicates the period of the fundamental frequency at 

the equilibrium; 

— on the right panel, the location of the stable equilibrium in the phase space (a, e, i) 
with the period in the color scale. 

The color code is the same for the left and right panels and it is truncated at the value 
of 100 years. For a larger period, we use the black color. 
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^ " ^3b''^j2 Semi-major axis [km] 



Fig. 9 Plot of the periods at the stable horizontal equilibrium (Eq. 11 with conditions 13), 
vertical equilibrium (Eq. 7 with condition 8) and (0, 0) equilibrium (with conditions 17) respec- 
tively in the upper, center and lower panels. The color code indicates the period (truncated to 
100 years) of the fundamental frequency at the equilibrium. On the left panels, the period with 
respect to the parameters (7, H^). On the right panels, the location of the stable equilibrium 
in the physical space (a, e, i) with its period. For the equilibrium (0, 0), e is always equal to 0. 
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6 Comparison of analytical and numerical solutions 

6.1 Comparison for all inclinations 

The analytical results of the simplified model described above are checked using a 
precise numerical integration of the complete set of equations of motion (1). For our 
test, we use Mercury's orbiter mission profile, which nominally puts the spacecraft into 
a high eccentric polar orbit. Numerical integrations were performed with the Bulirsch- 
Stoer (Stoer and Bulirsch 1980) integrator. We reproduce hereby afew characteristic 
plots of the numerical simulations to confirm our analytical theory (see Figure 10). 
Similar results have been obtained for a wide range of initial frozen orbit conditions. 

Figure 10 shows a very good agreement between analytical results and numerical 
simulations. 




Fig. 10 Comparison between analytical and numerical results. For the left and right panel, the 
initial conditions are ao = 6 407 km (7 = 1.000296), Oq = = 0°. In the left panel, for the 
lower orbit, we take eo = 0.545055 (G = 0.8384), iq = 76.646989° {H^ = 3.7492 X lO'^) and 
u)0 = 180°; for the upper orbit, we take eo = 0.6 (G = 0.8), jq = 78.221768° {H'^ = 0.26666666) 
and uiQ = 0° . For the right panel, the initial conditions are eo = 0.01 (G = 0.99995), «o = 
69.73104° (H^ = 0.1199999) and ljq = 0° . For the middle panel, the initial conditions are 
ao = 4 650 km (7 = 0.17096), eo = 0.3 (G = 0.9539392), iq = 77.89775° {H^ = 0.04) and 
Qq = Afo =1^0 =0°. The numerical model takes into account the contribution of J2 and 
C22 and the solar gravitational efEect, with starting epoch fixed at 14 September 2019. The 
analytical model is based on Equations (4, 5, 6). We plot the numerical integrations with 
continued lines and the analytical results with dashed lines. In the right panel, the numerical 
integration leads to a crash onto the planet. 



6.2 Comparison for polar inclination and explanation of the preliminary numerical 
results 

In the Figure 11, we present a graphical comparison between numerical integration and 
analytical results (contour plots of the Hamiltonian (3)) for an exact polar inclination. 
We see that the analytical theory is very close to the numerical integration for all initial 
eccentricities. We also notice that the addition of the C22 does not modify much the 
motion. 

In the right panel, we show two solutions close to the libration point and we see 
that, the closer the motion is to the libration equilibrium, the more the numerical 
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0.35 0.36 0.37 0.38 0.39 
k = e cos o> 

Fig. 11 Comparison between analytical and numerical results for exact polar orbiter. The 
initial conditions are ao = 6000 km (7 ~ 0.72), iq = 90° , J7o = 67.7°, u)o = -2°, M = 36.4°. 
The numerical and analytical model are the same of Figure 10. In dashed line the analytical 
result and in continued line the numerical integration. On the right a blow-up of the center of 
libration. 



integrations show a discrepancy with respect to the analytical results for the periherm 
libration: the frozen orbit of the analytical model shows no changes in eccentricity 
and argument of pericenter. On the contrary, the numerical orbit has short period 
oscillations but constant mean values of e and uj. 

Figure 11 allows us to explain the behaviors already seen in our preliminary numer- 
ical exploration (Fig. 1). In fact we can find there different orbits with a semi-major 
axis equal to 6 000 km corresponding to a vertical section in Figure 1. Then, on this 
section, we take some values of the eccentricity such that: 



— for e near to 0, in Fig. 1, we see a large value of the amphtude of variation of the 
eccentricity approximatively equal to 0.5 and a high value of the second derivative. 
In Fig. 11, for e equal to 0, we are on the separatrix. Therefore the eccentricity 
increases (roughly until 0.5) and a little shift of the initial eccentricity causes a high 
difference of the frequency. Thus the second derivative of the frequency is large; 

— for e close to 0.37, in Fig. 1, we see that the amplitude of variation of the eccentricity 
decreases until 0. 

In Fig. 11, ai e — 0.37, we find the stable point where the eccentricity is equal to 
a constant; 

— when e moves away from 0.37 to 0.5, in Fig. 1, we see that the amplitude of variation 
of the eccentricity increases from to 0.5 and for e — 0.5, the amplitude of variation 
of the eccentricity is maximal and the value of the second derivative is large. 

In Fig. 11, moving away from the equilibrium (e = 0.37) toward the separatrix 
(e ~ 0.5) we encounter larger and larger variations in e; 

— for e near to 0.58, in Fig. 1, we see that the amplitude of variation of the eccentricity 
is smaller than for e ~ 0.5. 

In Fig. 11, for e ~ 0.58, the pericenter circulates and the maximum of the amplitude 
of variation of the eccentricity is roughly equal to 0.58 — 0.12 = 0.46. 

— in Fig. 1, moving along the line e = 0, we pass from the region (F) to the region 
(E') at 5 577 km (Tab.l). For semi-major axis smaller than a = 5 577 km, we do 
not cross any separatrix and the amplitude of variation of the eccentricity is small. 
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Table 2 Comparison between the period of the equilibria determined in the analytical model 
and the period numerically obtained using NAFF. 



Initial condition 


Period [year] 


Error 


What 


a 


e 


I 


Analytical 


Numerical 


relative 


Equi. 


[km] 




[degree] 






% 


Kozai 


5 750 


0.4731 


58.37 


29.30 


29.25 


0.17 


Horiz. 


8 083 


0.4922 


77.68 


35.67 


35.61 


0.17 


Horiz. 


5 818 


0.5418 


71.93 


42.17 


42.26 


0.21 


(0,0) 


3429 


0.0 


47.64 


9.127 


9.135 


0.08 


(0,0) 


4 731 


0.0 


77.01 


56.594 


55.274 


2.38 



6.3 Frequency comparison 

To obtain a second independent validation of our analytical model, we numerically 
compute, using the NAFF algorithm (Laskar 1988, 2005), the period of the numerical 
solutions of the full system (1) obtained through numerical integration, and we compare 
it with the period of the equilibrium points of the simplified model. 

Table 2 provides a summary of this comparisons We can observe a very good agree- 
ment between the two methods. Some small differences can be explained as follows: 

— the exact equilibrium in the doubly averaged system is not the exact equilibrium 
in the full numerical model; 

— the full numerical model contains short period terms which disturb the long period 
dynamics. 



7 Discussions 

7.1 J2: the protector 

The aim of this section is to describe the protection mechanism of the coefficient J2 
on the increase of the eccentricity. We recall that our Hamiltonian (3), once we set the 
coefficient Ej^ = 0, reduces to the Kozai-Lidov Hamiltonian: 

In the Figure 12, we draw the possible phase spaces of this Hamiltonian. In the right 
panel {H'^ > 3/5) we have a similar behavior of our (E) case (Fig. 7). For the exact 
polar orbits {H^ = in the left panel of the Fig. 12), in the Kozai-Lidov Hamiltonian, 
all the probes are ejected: the eccentricity always grows up to 1. Instead, with the 
addition of the coefficient J2 we have the phase space (E') or (F) (Fig. 7) where it 
is possible that the eccentricity does not increase or that it remains at a fixed value. 
In the middle case {0 < H"^ < 3/5) we see that for an initial pericenter close to 0, 
the eccentricity increases. Instead, in our case, the phase spaces (A), (B), (C) and (E) 
(Fig. 7) show that it is possible to find initial condition (other than oj ~ ±7r/2) where 
the increasing of the eccentricity is naturally controlled. 



5(1-G0(1-|2 



. 2 
sin w • 



(21) 
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The J2 acts as a protection mechanism against the increase of the eccentricity due to 
the Kozai-Ltdov effect. This mechanism also appears for planets in tight binary systems 
(Saleh and F.A. 2009), where the general relativistic effects become dominant and can 
cause the periastron to precess on very short timescales. Therefore this precession can 
lead to the suppression of Kozai oscillations. 




270 270 270 

k=e cos (0 k=e cos ro k=e cos ro 

Fig. 12 All possible phase spaces for Kozai-Lidov Hamiltonian (21) with respect to the values 
of 



7.2 Local deformation of the Kozai-Lidov equilibrium 

We have seen that the condition to get the Kozai-Lidov equilibrium is (Eq. 8) 

57-1-5 

Actually there is a region where it is possible to find three real roots for G on a fonction 
of and 7. The conditions to have these three real roots are given by: 



864 0001/^^7^ + (^2 963 520//^^ - 1024iy^"^7*^ 



+ (1512 630//*^ - 13 965/?^ - 22 235 661/?^") 7^ -f 12 = 



and <^ 



We draw the solutions of this equation, denoted by KL^, that demarcates the region 
denoted (G), on the left panel of the Figure 13. Let us observe that this condition 
verified for large value of 7 (7 > 7203\/3/2) and for very small value of H'^ {H"^ < 
1/3087). An example of the phase space is plot in Figure 13 in the middle panels. In 
this region, the vertical Kozai-Lidov stable equilibrium bifurcates in two stable and 
one unstable vertical Kozai-Lidov equilibria producing thus a local deformation of the 
Kozai-Lidov equilibrium. We show an example of these three equilibria in the right 
panels of the Figure 13. Initial conditions close to these equilibria (external orbit in 
the right panels of Fig. 13) give rise to orbit librating around this set of three equilibria. 

It is possible to find that this bifurcation appears, in the (G) region, for a value 
of G smaller than y/3/21 ~ 0.082 478 6 corresponding to a value of the eccentricity 
e larger than \/438/21 ~ 0.996 59. Recalling the formula H — Gcosi, we obtain a 
minimal inclination of 87.27°. 
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Fig. 13 Local deformation of the Kozail-Lidov equilibrium. The bifurcation lines in the left 
panel with the new region (G) demarcated by the two curves KL3. Example of generic contour 
(for (G) region) of the Hamiltonian (3) in (fc, h, i) space in the middle panels. A zoom of the 
local deformation in the right panels. 



7.3 J3 discussion 

In Paskowitz and Scheeres (2006) the authors included the J3 (the "pear shape" of 
the central body) Europa's effect in their system. They noticed that the coefficient J3 
caused an asymmetry between the solutions of the frozen orbits for uj — ±7r/2 but they 
did not explain the reasons of this beavior. 

The potential arising from a central body with a J3 7^ is given by 

The averaged Hamiltonian is then 

SGMJsRf, . . 5 



2a4(l-e2)5/2 



e sin u} sin til sin i 



Using our variables G — \/l — e^, H = Gcosj, we can define the dimensionless (divided 
by GM/a) potential that we can add to the Hamiltonian (3): 



JsRp 



yi-G2 sino; - i/2 (5//2 _ ^: 



2\ 
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Introducing the coefficient 6 = = , the equations of motion (6 and 5) can be 
rewritten in compact form as follows: 

(G = Fi(G, J?,7) sinwcosw + F2{G,H,S) cosw 

\w =F3(G,i?,7) + F4(G,-ff,7) sin^w + F5{G,H,5) sinw 

where the functions Fi, F3 and F4 can be easily identified in equations (5) and (6). 
The functions F2 and F5 come from the J3 effect and they are proportional to 5. 

7.3.1 Vertical equilibria - Kozai-Lidov equilibria: cosw = <^ a; = ±7r/2 

Let us observe that the addition of J3 effect causes an asymmetry in the frozen orbit 
solutions not present before. Indeed, for uj = -k/2 the condition of equilibrium is given 

by 

FsiG, H,^) + FiiG, H, 7) + F5(G, H,S) = 
whereas for oj = —it/2 the condition of equilibrium is given by 

F3{G,H,'y) + Fi{G,H,^) - F5{G,H,S) = 0. 

Then, for a small coefficient S, the asymmetry is not important. However for a large 
value of this coefficient, the asymmetry could be important until the elimination of one 
of two equilibria. 

7.3.2 Horizontal equilibria. 

For horizontal equilibria, the condition of equilibrium (G = 0) becomes: 

Fi (G,i?, 7) sin w + F2 (G, i?, <5) = <(=^ sinw = -F2/F1 -e. 

Then the "horizontal" equilibria appear for non-zero values of the pericenter cj = —e 
and w = TT + e. The condition to obtain w = becomes: 

FsiG, H, 7) + F4(G, H, 7)£2 - F^iG, H, S)e = 0, 

that induces a shift in the equilibrium in G and u) variables with respect to the case 
"J2 + third body". 

7.3.3 Modifications of the phase space 

For illustration, in Figure 14, we draw the contour plots of the new Hamiltonian for 
different values of J3 (or for different values of S). We see that when the 6 coefficient 
increases (in absolute value), the vertical equilibrium goes down while the horizontal 
equilibrium goes below the "line smuj = 0". We point out that, from some values of S, 
the equilibrium co = — 7r/2 disappears (Fig. 14 right panel). 
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Fig. 14 Distortion of the phase space (for a Mercury's orbiter) due to J3 effect. The initial 
conditions are a = 5 900 km (7 ~ 0.66) and = 0.06. J3 (respectively 5) is equal to (0), 
-J2/IO (-0.041356) and - J2/2 (-0.20678) in left, center and right panels. 



7.3.4 BepiColombo and other missions 

At present time, the semi-major axis of the two orbiters (MPO & MMO) of the Bepi- 
Colombo mission are respectively equal to 3 394 km and 8 552 km. The MPO altitude 
corresponds to our (E') phase space where the eccentricity vector has a circular concen- 
tric motion. The MMO initial conditions, without thrust correction, leads to a crash 
onto the Mercury surface after 3 years. Thanks to our theory, we can choose another 
initial condition a = 7355 km and e = 0.652, that avoids the crash on Mercury and 
whose eccentricity vector is fixed. 



8 Conclusions 

The orbit dynamics of a space probe orbiting a planet or a natural planetary satellite 
has been investigated. The proposed model includes the effects of J2 for the central 
body and the perturbation of the third body. We have developed a doubly averaged 
Hamiltonian and studied the location of the stable and unstable frozen orbits. Our 
analytical approach allows us to compute also the periods of the free librations at 
the equilibria. The analytical results have been checked and validated numerically 
by performing numerical integrations of the complete systems. Our theory is able to 
explain the behavior of our preliminary numerical investigations where the variation 
of the amplitude of the eccentricity is null and the presence of a separatrix has been 
found by numerical investigation. The theory is general enough to be applied to a wide 
range of probes around any planet or any natural planetary satellite, provided that 
they respect the hypotheses used to obtain our Hamiltonian model. 

We have shown the protection mechanism of the coefficient J2 on the increasing 
of the eccentricity due to Kozai-Lidov effect. This mechanism is therefore able to find 
a larger number of frozen orbits than for the only Kozai-Lidov problem. We have 
also explained the asymmetry of the frozen equilibria caused by the addition of the 
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coefficient J3. Wc have also brought to the hght a local deformation of the Kozai-Lidov 
equilibria that appears at high eccentricity, high inclination and large value of 7. 

It would be interesting to take this theory into account to choose the intial semi- 
major axis and eccentricity of an orbiter for future missions around planets or planetray 
satellites. 
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Appendix: "Le theoreme d'algebre de Sturm" 

Let f(x) be a polynomial of positive degree with real coefficients and let {fo{x), fl {x), fiix), . . . , /s(a;)} 
be the standard sequence for f{x) such as 





= nx) 






Mx) 


= fix) 






h{x) 


= <lo{x)fi{x) - 


Mx), 


deg /2 < deg /i 


h{x) 


= H(x)f2ix) - 


hix), 


deg fs < deg /2 


h+\{x) 


= qi-i{x)fi{x) 


- fi-i{x) 


, deg /i+i < deg ft 


until /s+i(a;) 


= 







where fi—i is obtained by the Euclidean division: /i+i = qi—ifi — fi—i- Assume [a,b] is an 
interval such that /(a) ^ ^ f{b)- Then the number of distinct roots of f{x) in [a,b] is 
Va — VJ, where Vc denotes the number of variations in sign of {/o(c), /i(c), . . . , /s(c)}. The 
are dropped from the sequence. 



